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I. INTRODUCTION 

Modern probes of material properties, such as the 
new inelastic neutron scattering facilities, are reaching 
such unprecedented sensitivity that they can measure the 
spectrum not only of a single quasiparticle excitation, 
but even two-particle excitations These quasiparti- 
cles can collide, scatter, or form bound states just like 
elementary particles in free space. The spectrum of the 
multiparticle excitations is a crucial indicator of the un- 
derlying dynamics of the system. 

One of the principal theoretical means of predicting the 
excitation spectrum is the method of high-order pertur- 
bation series expansions (2|. We have previously used a 
'linked-cluster' approach to generate series expansions for 
2-particle states in 1-dimensional models @, but for 2- 
dimensional models the only high-order calculations car- 
ried out so far have been those of Uhrig's group (e.g. [H), 
using the 'continuous unitary transformation' (CUTS) 
method, which is of only limited applicability. One of 
our aims here is to extend the linked-cluster approach to 
2-dimensional models, starting with the bilayer model as 
a simple example. 

The 5 = 1/2 bilayer Heisenberg antiferromagnet has 
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FIG. 1: The bilayer Heisenberg model on a square lattice. 



attracted continuing interest from both experimentalists 
and theoreticians. Experimentally, it is of interest be- 
cause many of the cuprate superconductors contain pairs 
of weakly coupled copper oxide layers @, H, S @|- Re- 
cently, the organic material piperazinium hexachlorod- 
icuprate has also been found to have a bilayer structure 
Q . Theoretically, it is of particular interest because it is 
one of the simplest two-dimensional systems to display 
a dimerized, valence-bond-solid ground state, when the 
interplane coupling is large. There have also been dis- 
cussions of the model in the presence of a magnetic field 
0, or doping @, El El, or disorder [H|. 

The structure of the model is shown in Figure [TJ with 
5 = 1/2 spins on the sites of the lattice, and Heisenberg 
antiferromagnetic couplings Ji between the planes, J\ 
within each plane: 

l = 1,2 <i,j> i 

where I = 1,2 labels the two planes of the bilayer. The 
physics of the system then depends on the coupling ratio 
A = J1/J2. At A = 0, the ground state consists simply of 
5 = dimers on each bond between the two layers, and 
excitations are composed of 5 = 1 'triplon' states on one 
or more bonds. At large A, where the J\ interaction is 
dominant, the ground state will be a standard Neel state, 
with 5 = 1 'magnon' excitations. At some intermediate 
critical value A c , a phase transition will occur between 
these two phases. It is believed that this transition is 
of second order, and is accompanied by a Bose-Einstein 
condensation of triplons/magnons in the ground state. 

Theorists have discussed this model using series expan- 
sion methods [H], EEEH> quantum Monte Carlo )QMC) 
simulations at small temperatures [IE El) EH , Schwinger- 
boson mean-field theory [2(| SII, and spin- wave theory 
[H SI Sip!. The QMC analysis of Sandvik and 
Scalapino [19| found the transition at A c = 0.398(3), with 
a critical index v]simeq§.7. in agreement with the 0(3) 
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nonlinear sigma model prediction, while the exponent- 
biased series analysis of Zheng [l| put the critical point 
at A c = 0.394(1). Early spin-wave estimates ;23] were 
well away from this position, but the improved Brueck- 
ner approach of Sushkov et al. [24, 25] gave a remarkably 
accurate estimate of the critical point and critical index, 
and also the 1-particle dispersion in the model. 

Our particular aim here is to study the two-triplon 
states within the dimerized regime, with particular em- 
phasis on the occurrence of bound states, and to explore 
their behaviour in the vicinity of the critical point. The 
two-particle bound states can give important insights 
into the dynamical behaviour of the model. It is also 
possible that they may be detected experimentally at the 
new generation of inelastic neutron scattering facilities, 
or by other means. 

We use two methods to investigate the two-particle 
states. A modified triplet-wave approach, described in 
Section HH gives a qualitative picture of these states, 
valid at small couplings A. Series expansion calculations, 
sketched in Section IIII) are then used to obtain more 
accurate results, and to explore the behaviour near the 
critical point. Series expansions are also presented for 
the single-particle and total transverse structure factors. 
Our conclusions are summarized in Section HVl 

II. MODIFIED TRIPLET- WAVE THEORY. 

Analogues of spin-wave theory in a dimerized phase 
have been discussed by several authors. Sachdev and 
Bhatt [26| used a 'bond-operator' representation to de- 
scribe the dimers and their spin-triplet excitations, which 
employed both triplet and singlet operators, with a con- 
straint between them to ensure that no two triplets can 
occupy the same site. The constraint is awkward to im- 
plement, and so Kotov et al. [25[ discarded the singlet 
operator, and replaced it by an infinite on-site repulsion 
between triplets, implemented via a self-consistent Born 
approximation, valid when the density of triplets is low. 
We have presented an alternative approach [271 ] . where 
the exclusion constraint is implemented automatically by 
means of projection operators. The absence of any con- 
straint makes the formalism easier and more transparent 
to apply, but at the price of extra many-body interaction 
terms. This is the method used here. 

The Hamiltonian for the Heisenberg bilayer systemn 
can be rewritten 

h = s ii • s 2i + a £ Sa; ■ Sa j (2) 

i 1=1,2 <i,j> 

For A = 0, the system reduces to independent dimers as 
shown in Figure [TJ Let us consider a single dimer with 
two spins Si,S2- The four states in the Hilbert space 



consist of a singlet and three triplet states with total 
spin S = 0, 1 respectively, and eigenvalues 

S i' S 2- { +1/4 (S = l) (3) 
We denote the singlet ground state as |0), and introduce 
triplet creation operators that create the triplet states 
out of the vacuum |0), as follows 
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Then the spin operators Si and S2 can be represented 
in terms of triplet operators by 

S± a = 2^ck(1 ttftryj + (1 tl f ty)t a teapot pt*y\ 

S 2a = -[-4(1-^)- (l-t+* 7 )* Q 

itaft-ftjjt^] (5) 

where a, /?, 7 take the values x, y, z and repeated indices 
are summed over. This is similar to the representation 
of Sachdev and Bhatt [26[ , except that we have omitted 
singlet operators s', s, but used projection operators (1 — 
tiyt-y) instead. Assume the triplet operators obey bosonic 
commutation relations 

[t a ,t'p] = S a p, (6) 

then one can show that within the physical subspace (i.e. 
total number of triplet states is or 1), the representation 
([5|) obeys the correct spin operator algebra 

[Sia, Sip] = iCaPfSiy, [52a; £2/3] = i£aPfS2y, (7) 



[Sl a ,S 2 p]=0 (8) 

S? = S 2 2 = 3/4, S x • S 2 = *tt Q - 3/4 (9) 

The projection operators ensure that we remain within 
the subspace. 

Returning to the bilayer system, we can now define 
triplet operators t' na , t na for each dimer n in the system. 
For a system of N dimers, the Hamiltonian now can be 
expressed in terms of triplet operators as 
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3iV A A A A A A A r / f t | \ / t 

" = ^ I" 2^ t na tnu + ^ 2^ {t ia tj a + t ia tj a + ti a tj a + t ia tj a \ — — {(^cA,3^/3 + tiptif3tia){tj a + Aj 

n <y> <ij> 

A 



+ ^ {(4a 4/3**/3 + 4/3^/3^ Q ) (447^7 ^T^T^ ^ 

I 



(10) 



<ij> 



This expression includes terms containing up to 6 triplet 
operators. 

Next, perform a Fourier transform 



4 = (^) 1/2 E^ k ' n 4 



(ii) 



t ka = (^) 1/2 E e4k - n ^ 



(we set the spacing between dimers d = 1), then the 
Hamiltonian becomes 



3iV - — 

H = +Y i tl a t ka + \Y I 7k[*L* t - kQ +tkai-ka + 24 Q i kQ ] 

k k 

-4 X)^ 1+2 + 3 - 4 [^l«*2«4 7 *47 + 4 7 *3 7 *2a<la)(7l + 72)] + 5l+2-3-4[*L4 7 *3a*4 7 (7l + 72 + 73 + 74) 



+7l-3i 1( 3i 2 /3^3 7 *4 7 — 7l-4^ 1/ 3*2 7 ^3/3i4 7 )]} 

+ ~j\f2 X!'f <5l + 2 + 3 + 4 - 5 - 6 [ 73 + 4 - 6 ^la474a4 / 3 < 5 7 i6/3 + 4/34 7 *4/3*3a*2 7 *la)] 
1-6 

+^l+2+3-4-5-6il a *2/3*3 7 *4 7 i5/3i6a(73-4-6 + 72+3-4)} 



(12) 



where the indices 1 — 6 are shorthand for momenta ki 
k6, and 



7k = ~{cosk x + cos k y ) 



(13) 



serves the boson commutation relations 



(15) 



for the square lattice. Henceforward, we drop the 6- and is intended to diagonalize the Hamiltonian up to 

particle terms quadratic terms. After normal ordering, the transformed 

Finally, as in a standard spin- wave analysis, we per- Hamiltonian up to fourth order terms reads 
form a Bogoliubov transform 



tlaa = CkTka + S^T^^ (14) 

where Ck = cosh6>k, Sk = sinh6>k, 0-k = $k, which pre- where the constant term is 
I 

Wo = 3iV | ~ + R 2 + 2\(R 3 + R A ) - 2\[2{R 3 + i? 4 )(i?i + 4i? 2 ) + ^ 7i-2(cisic 2 s 2 - s 2 s 2 )] 

I 12 

+2A[(i? 3 + R 4 )(Ri + m 2 f + ^3 E 71+2-3 (ciSi(4c 2 S2C 3 S3 + 6c 2 s 2 s 2 + 6s 2 s 2 ) + 4s 2 s 2 s 2 )] 



H = W Q + H 2 + H 3 + H 4 , 



(16) 



(17) 
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expressed in terms of the momentum sums 

Ri = jjJ2 CkSk 



4 



N 4^ 



CkSk7k 



where 



(18) 



The quadratic terms are 

H 2 = ^[£ k T kQ r kQ + Qk(r ka r_ kQ + r ka r! kQ )] (19) 



£k = 



(4 + «k)(l + 2A 7k) + 4A 7k c k Sk 

-A[4(c k + 4)(7k(i?i + 4i? 2 ) + 4(i? 3 + i?4)-^E «?Tk-i) 



+8c k s k ( 7k (i? 1 + 4i? 2 ) + (i? 3 + i? 4 ) + jz 2 cisi7k-i)] 



(20) 



Ok = 



c k s k (l + 2A 7k ) + A 7k (4 + 4) 

-A[2(4 + 4)(7k(i?i + 4R 2 ) + (i? 3 + Ik) + ^ 51 c i s i7k-i) 



+4c k s k ( 7 k(i?i + 4i? 2 ) + 4(i? 3 + i? 4 ) - 53 s i7k-i 



The fourth-order terms are 



(21) 



1+2+3+4 



$ 4 (n t a T 2a' r 3 t 7 7 47 + r l Q T 2a T 37 T 47 ) + 5l+2-3-4($4^ T la r 2a r 3 7 T 4 7 + ^ Wa^V 3 " 7 ^) 



( 2 U _t 



s( 3 U J 



1234 



+ '5l+2+3-4$i 4) (Tl t Q ' 7 '2 t c e T 3 t 7 T4 7 + T 4 V 3 f r 2a T la)] 



where we have used the shorthand notation 1 • • • 4 for mo- 

(i) 

menta k\ ■ ■ ■ fc 4 , and the vertex functions $> 4 arc listed 
in Appendix A. These results were obtained or con- 
firmed using a symbolic manipulation program written 
in PERL. 

The condition that the off-diagonal quadratic terms 
vanish is 

Qk = 0. (23) 
In a conventional spin-wave approach, this would be im- 



(22) 



I 

plcmcntcd in leading order only, giving the condition 
2s k c k 2A 7k 



tanh 2# k 



[1 + 2A 7k ) 



(24) 



This would leave some residual off-diagonal quadratic 
terms, arising from the normal-ordering of quartic op- 
erators. In a 'modified' approach (28|, we demand that 
these terms vanish entirely up to the order calculated, 
giving the modified condition 



J 



tanh20 k = — 



2A[ 7k - 2(7k(fli + 4fi 2 ) + (R 3 + R^) + w Si ciaiTk-i)] 
[1 + 2A( 7k - 2( 7k ( J R 1 + 4i? 2 ) + 4(i? 3 + Ri) ~ £ Ei 47k-i))] 



(25) 



Self-consistent solutions for the N equations (|25p . with easily be found by numerical means, starting from the 
the four parameters R\ ■ ■ ■ i? 4 given by equation (fig)) , can conventional result (EM)) . 
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FIG. 2: Perturbation diagrams contributing to the ground- 
state energy. 



A. Expansion in powers of A 

As a first check on the formalism, one may calculate the 
leading terms in an expansion of the energy eigenvalues 
in powers of A. Solving the modified equation (|25|) self- 
consistently to order A 2 , we find 



A 2 

«k = -A7k + y( 4 7k 
c k = 1 + ^A 2 7 £ 



7k -1) 



(26) 



with the lattice sums ([T8|> 



Ri 



0{\ 4 



i?2 = -r 



A 2 
4 



y+0(A 4 ) 



A A 2 
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+ 0(A 3 ), i? 4 = 0(A 3 ) (27) 



The leading-order behaviour of the vertex functions may 
easily be deduced from Appendix A. 

Substituting in equation (|17[) . the ground state energy 
per dimer is 



£n = ~ —3 A 

N L 4 4 8 J 







(28) 



in agreement with dimer series expansion results previ- 
ously obtained for this model [HI]. One can easily show 
that perturbation diagrams such as those in Figure [5] do 
not contribute until 0(A 4 ) or higher. 

The energy gap at leading order can be found from 
equation (f2"D|) : 



E k ~ 1 + 2A 7k + 4A 2 - 2A 2 7 £ A 







(29) 



Note that in linear spin-wave theory, when tanh 2# k is 
given by (|24[) and the energy gap is given by the first line 
of equation (|2H)) . the energy gap is 



yfi + 4A 7k 



(30) 



which vanishes at A = 1/4, 7k = — 1) i.e. k = (-7T, ir). This 
marks a phase transition with critical index v — 1/2, and 
the end of the dimerized phase, in this approximation. 





FIG. 3: Perturbation diagrams contributing to the one- 
particle energy. 



The perturbation diagram Figure [3^) also contributes 
to the energy gap at order A 2 . Note that diagram[3|i)does 
not appear in the formalism of Shevchenko et al. |24l . [25[ ; 
the extra terms in our formalism are needed to implement 
the hardcore constraint that two triplons cannot occupy 
the same site. At leading order, the contribution of this 
diagram is 



AE] 



3a) 



-2A 2 



A^O 



(31) 



(see the next section for further details). This gives a 
total single-particle energy 



l + 2A 7k + 2A 2 (l- 7 2 ), A^O (32) 



which again agrees with series expansion results [jq ]. 

The minimum energy gap lies at k = (7T,7t). If we 
compare equation (]32[) at small momentum p = (tt, n) — k 
with the continuum dispersion relation for a free boson, 



€k ~ y m 2 c 4 + p 2 c 



(33) 



we readily discover the leading behaviour of the effective 
triplon parameters, i.e. the triplon mass 



m~ -[1 -2A + 0(A 2 )] 
A 

and the 'speed of light' or triplon velocity 



A + 0(A 3 ) 



(34) 



(35) 



in lattice units. Note that the mass diverges and the 
speed of light vanishes as A — > 0. 
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-0.75 



-0.76 



«> -0.77 




and 



-0.78 - 



-0.79 - 



V = Ha 



(38) 



(6-particle terms being neglected) we can treat Hq as 
the unperturbed Hamiltonian and V as a perturbation 
to obtain the leading-order corrections to the predictions 
for physical quantities outlined in the previous section. 
Numerical results for the model have been obtained us- 
ing the finite-lattice method. The momentum sums are 
carried out for a fixed lattice size L x L — N, using cor- 
responding discrete values for the momentum k, e.g. 



0.25 



FIG. 4: Ground-state energy per dimer as a function of A. 
The solid line is the estimate from series expansions, and the 
dashed line is the triplet-wave estimate. 



n = 1 , • • • L 



2im 

x = — , 

27rm 

y = m = 1, • • • L 



(39) 



B. Numerical Results 



Writing the Hamiltonian as 



Results were obtained for L up to 100. 



where 



H = H n 



V 



H = Wo + H 2 



(36) 1. Ground-state energy 

The leading correction to the ground-state energy cor- 
responds to the diagram in Figure [2^,). Its contribution 

(37) is 



Ac 



2a) 



-3A 2 v 

N 3 1^ 

1234 



J l+2+3+4 



$i 1} (1234) 



[E x +E 2 + E 3 + E 4 ) 



3$^ (1234) + $^(1324) + $^(1423) 



(40) 



r 



In leading order one can show that this term is 0(A 4 ), 
whereas diagrams such as Figure[3b) are 0(A 6 ) or higher. 
Figure U shows the behaviour of the ground-state energy 
as a function of A resulting from this modified triplon 
theory, as compared with the high-order dimer series cal- 
culations of Zheng [HI]. It can be seen that out to A ~ 0.1 
there is quantitative agreement between our calculation 
and the series estimates, but some discrepancy emerges 
at larger A. 



2. One-particle spectrum 



The leading correction to the one-particle spectrum 
corresponds to the diagram in Figure [3^,) . Its contribu- 
tion is 



AE] 



3a) 



1+2+3-k 



123 



$i 4) (123/fc) 



(Ek — Ei — E 2 — E 3 ) 



r 



3$i 4) (123/fc) + $r (321*0 + $r(312fc) 



i(4) 



(4), 



(41) 



In leading order, this term is 0(A 2 ), as stated in the higher, 
previous section, while diagrams like [3b) are 0(A 4 ) or 
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(0,0) (jt.O) (tt.tt) 



(0,0) 



FIG. 5: [Colour online] One-particle dispersion relation at the 
critical point (y — 1/A), as estimated from both dimer (solid 
line) and Ising (dashed line) expansions [151 ]. 



1.5 



Triplet wave, 2p terms — i — 
4p terms 
4p with corrections * 

Brueckner approach — s — 

Series expansion — ■ — 




FIG. 6: Energy gap at k = (n, n) as a function of A. The solid 
squares show the series estimates I15I1 , and the open squares 
are results from Shevchenko et al. [24j, while the stars show 
the improved triplet-wave results. The contributions from 2- 
triplon and 4-triplon terms are shown separately. 



the triplet-wave estimates begin to diverge, as higher- 
order contributions become more important. The self- 
consistent Born approach of Kotov et al. [24], [H| is more 
accurate than our approach at large A; but neither ap- 
proach can compete with scries methods for accuracy. 
Our object here mainly is to understand the qualitative 
behaviour of the model. 



3. Two-triplon bound states 

It has been found in previous studies of dimerized an- 
tiferromagnetic systems in one dimension 0, [29[ that 
the quartic terms in the Hamiltonian lead to attraction 
between two elementary triplons, giving rise to 5 = 
and 5=1 bound states. We look for solutions of the 
two-body Schrodinger equation 



H\^) = E\i>). 



(42) 



The two-body wave functions |i/>(K)) can be written 
as follows: 

Singlet sector (S = 0): 

where K is the centre-of-mass momentum and q the rela- 
tive momentum of the two particles, and the scalar wave 
function is symmetric, 



^ S (K,-q) = ^ S (K,q) 
Triplet sector (S = 1 ): 



(44) 



IV£(K)) 4£ ^ 7 ^ T (K,q)4 /2+q/3 4 /2 _ q7 |0) 

(45) 



with the wave function antisymmetric 

V> T (K,-q) = -^ T (K,q). (46) 

We will not write out the quintuplet states explicitly. 

From equation (j42f one can readily derive the inte- 
gral Bethe-Salpeter equation satisfied by the bound-state 
wave functions: 



The dispersion of the one-particle energy as a function 
of momentum k at the critical point is illustrated in Fig- 
ure [5j as estimated from two different series expansions 
by Zheng [l5j. It can be seen that the two expansions 
agree well at the critical point, and that the energy gap 
vanishes there at the Neel point k = (tt, tt). 

The triplet-wave and series estimates of the energy gap 
at k — (tt, tt) are compared in Figure [S] It can be seen 
that the inclusion of the corrections from diagram^) im- 
proves the agreement with series substantially, bringing 
quantitative agreement out to A ~ 0.15. Beyond that, 



[E S ' T 'Q(K) - £ K/2+q - £ K / 2 - q ]V^ T ' Q (K,q) = 

l^M^(K,q,p)^ T ' Q (K,p) (47) 
p 

in each sector S,T or Q. 

In leading order, the scattering amplitudes 
M ,T,< ^(K,q, p) are simply given by the 4-particle 
vertex from the perturbation operator V, Figure [7^,). 
Hence we find for the different sectors: 



M H (K, q, p) = -2A[3$ 4 2) (K/2 + q, K/2 - q, K/2 + p, K/2 - p) + $ 4 JJ+ (K/2 + q, K/2 - q, K/2 + p, K/2 - p)] (48) 



(3)+/ 



M T (K, q, p) = -2A$i 3) "(K/2 + q, K/2 - q, K/2 + p, K/2 - p) 



(49) 



M «(K, q, p) = -2A$f )+ (K/2 + q, K/2 - q, K/2 + p, K/2 - p) 



(50) 




+ P 



^4 , ^4 



(a) 



; (4) 



X 




^4 



t 



(c) 



X^4 (2) ^4 <3> 





(d) 



FIG. 7: Perturbation diagrams contributing to the 2-particle 
scattering amplitude. 



where the wave function is once again symmetric in the 
quintuplet sector 



V> Q (K,-q)=V> Q (K,q) 



(51) 



and the symmetric and antisymmetric pieces of the ver- 



fS) 

tex function <I> 4 are defined: 



(3)± _ 1 [fh (3) 



[$^(1234) ±$ 4 3) (1243)]. (52) 



At leading order in A, we find 

Af 5 (K,q,p) 2A[ 7p+q + 7p _ q + 7K/2+P 

+7K/2+q + 7K/2-p + 7K/2-q] (53) 




(0,11) 



(n,Jt) (2jt/3,2it/3) 



FIG. 8: [Color online] Dispersion relations for the two-particle 
bound/ ant ibound states at A = 0.1, along symmetry lines in 
the Brillouin zone, as calculated from the triplet-wave expan- 
sion. The 2-particle continuum region is shaded. 



M T (K,q,p)~ A[ 7q+ p- 7q -p] 



(54) 



and 



Af Q (K,q,p) ~ A[7 q+ p+7 q _p-2( 7 K/2+p 

+7K/2+q + 7K/2-p + 7K/2^q)] (55) 

Then restricting ourselves to the particular momen- 
tum K = (tt, 7t), simple solutions to the Bethe-Salpeter 
equation (|4"T|) can be found: 



* S ' Q (K,q) - (cosq x ±cosq y ) 
* T (K,q) ~ (smq x ±smq v ) 



(56) 



corresponding to nearest-neighbour pairs of triplon exci- 
tations, with energies: 

E S {K) - 2- A 

E T {K) ~ 2 -A/2 

E Q (K) ~ 2 + A/2 (57) 



9 



2 



in 



0.5 





f 


l 


\ 

. \ 




_ jyj^o.i 

rJ/J^O.3942 


J,/J 2 =0.1, 


.2, 0.3, 0.3942 



(0,0) (71,0) (tt.tt) (0,0) 



FIG. 9: [Color online] The total static structure factor S(k) 
as a function of k, for various couplings A = 3\ / J2 . 

Since the 2-particle continuum is confined strictly to 
E C ont = 2 at this order and this momentum, we see that 
the singlet and triplet states are bound states lying below 
the continuum, while the quintuplet states are antibound 
states lying above the continuum. There are two degen- 
erate states in each case, corresponding to the ± signs in 
equation (fijrj]) . or to the two possible axes x and y of the 
nearest- neighbour pairs. At higher orders these states 
may mix and separate. 

These results are easily understood in a qualitative 
fashion. For an S z — 2 excitation, for example, the spins 
on the nearest-neighbour sites are necessarily aligned par- 
allel, giving rise to a repulsive interaction; whereas for 
S — or 1 the neighbouring spins can be aligned either 
parallel or antiparallel, allowing the possibility of an at- 
tractive interaction. 

Solving the wave equation (|47[) with vertex functions 
given by the leading order approximations (|48[) - ([50]) . we 
obtain numerical solutions for the 2-particle spectrum, as 
illustrated in Figure [H at a coupling A = 0.1, near mo- 
mentum k = (ir, tt). It can be seen that the pairs of de- 
generate S = and S — 2 bound/antibound states split 
as one moves away from (jr, tt), and all states eventually 
merge into the continuum. 

III. SERIES EXPANSIONS 

We have performed a standard dimer series expansion 
@; H3| for this model, where the Hamiltonian is written 



as 

H = H a + XV (58) 
H = ^S li -S 2i (59) 

i 

v = E E s » • s u ( 6 °) 

1=1,2 <ij> 

and perturbation series are generated for the quantities 
of interest in powers of A, using linked cluster methods. 
Details of the linked cluster approach are reviewed in 
• In very brief summary, the ground-state energy per 
dimer can be written as a sum of the irreducible contri- 
butions (cumulants) coming from every connected clus- 
ter of dimers which can be embedded on the lattice, the 
order of the contributions rising with the size of the clus- 
ter. The 1-particle energies can be written in terms of 
irreducible transition amplitudes A\(i,j) of the effective 
Hamiltonian [1(1 , which consist of a sum over all linked 
clusters connected to i and j, the initial and final po- 
sitions of the 1-particle excitations. And finally, the 2- 
particle energies can be written in terms of irreducible 
transition amplitudes A2(i, j; k, I) of the 2-particle effec- 
tive Hamiltonian consisting of a sum over all linked 
clusters connected to (i,j) and (k, I), the initial and final 
positions of the 2-particle excitations. The amplitudes 
A2 are then employed in the 2-particle Schrodinger or 
Bcthe-Salpctcr equation to calculate the energy for as a 
function of momentum. We use a finite-lattice approach 
[2| for this purpose, where the Schrodinger equation is 
solved on a finite lattice in position space, of sufficient 
size to ensure convergence of the results. 

Once a perturbation series in A has been calculated 
for a given quantity, it can be extrapolated to finite A 
using Pade approximants or integrated differential ap- 
proximants. 

Zheng [l5[ has previously calculated series for the 
ground-state energy and 1-particle excitations. These re- 
sults have already been compared with the triplet-wave 
predictions in Figures 2J [5] and [5] 



A. Structure Factors 

Figures 191 and ITD1 show some series results for structure 
factors, which have not been shown before. Figure [9] 
shows the total static transverse structure factor S(k) = 
S~* (k) as a function of k at various couplings A = J1/J2, 
where (k) is the Fourier transform of the correlation 
function: 

5+ ~( k ) = jf E e * k ' (ri ~ rj) < s l s i >° ( 61 ) 
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TABLE I: Series coefficients of A in the expansions for the 1-particle structure factor Si p and integrated structure factor S 
at momenta k = (ir,n) and (0,0). 



N 


Si p (tt,ty) 


S(n,n) 


5i p (0,0) 


5(0,0) 





1 .00000000000000D+00 


1.00000000000000D+00 


1 .00000000000000D+00 


1.00000000000000D+00 


1 


2.00000000000000D+00 


2.00000000000000D+00 


-2.00000000000000D+00 


-2.00000000000000D+00 


2 


5.00000000000000D+00 


5.43750000000000D+00 


3.00000000000000D+00 


3.43750000000000D+00 


3 


1.20000000000000D+01 


1.24375000000000D+01 


-7.00000000000000D+00 


-6.56250000000000D+00 


4 


2.60000000000000D+01 


2.73476562500000D+01 


1.42500000000000D+01 


1.48476562500000D+01 


5 


6.19609375000000D+01 


6.16328125000000D+01 


-3.08359375000000D+01 


-3.09609375000000D+01 


6 


1.45859863281250D+02 


1.46245605468750D+02 


6.65551757812500D+01 


6.68159179687500D+01 


7 


3.60063964843752D+02 


3.57834899902344D+02 


-1.51234863281252D+02 


-1.51278381347656D+02 


8 


8.71365653991730D+02 


8.80394332885743D+02 


3.23292167663603D+02 


3.28300582885742D+02 


9 


2.13146787007666D+03 


2.15030324554441D+03 


-7.25282606760795D+02 


-7.27275304158507D+02 



All results are for k z = tt, probing intermediate states 
antisymmetric between the planes, and we only refer to 
k = (k x ,k y ) hereafter. 

The dominant feature is a large peak at the Neel point 
k = (tt, tt), which appears to become divergent as A — > A c . 
This behaviour is qualitatively very similar to that seen in 
the alternating Heisenberg chain (AHC) in one dimension 
[3l1 |. Figure [10] shows the ratio of the 1-particle structure 
factor Sip(k) to the total S(k) as a function of k. The 
1-particle contribution generally remains the dominant 
part of the total, particularly near the Neel point. This 
behaviour is again reminiscent of the AHC [31| ■ 

Further information may be obtained from the series 
for S(k) and Si p (k) at the Neel momentum (jr, tt), which 
are given in Table I. A Dlog Pade analysis of these series, 
biased at A c = 0.3942, shows both 5(k) and Si p (k) di- 
verging as A — > A c with exponents —0.68(1) and —0.69(1) 
respectively. The series for the ratio Si p /S shows no sign 
of a singularity at this point, remaining almost constant, 
within 2% of unity at all couplings. This behaviour is 
quite different from the AHC case [35(, where the ratio 
vanishes logarithmically at the critical point. 

These results should be compared with theoretical ex- 
pectations. From scaling theory (see Appendix B), the 
1-particle structure factor in the vicinity of the critical 
point Si p (n, tt) should scale like (A c — A)^ -1 ^, at the 
critical (Neel) momentum. For the total structure fac- 
tor at this point, scaling theory gives exactly the same 
exponent (see Appendix B) . We expect this transition to 
belong to the universality class of the 0(3) model in 3 di- 
mensions, which has critical exponents [32| v = 0.707(4), 
r] = 0.036(3), hence we expect (rj — l)v = —0.682(5), 
which is quite compatible with the numerical estimates 
obtained above. 

How does Sip behave at the critical coupling away from 
the Neel momentum? In the transverse Ising model [33j , 
it was found that the 1-particle residue function A(k) (see 



CO 
a 

U2 




(0,0) 



FIG. 10: [Color online] The ratio 5i p (k)/S(k) of the 1-particle 
static structure factor to the total static structure factor as a 
function of k, for various couplings A = J1/J2. 



Appendix B) vanishes like (A c — X)' r,L ' at all momenta, with 
a small exponent rjv — +0.025(3), so that Si p vanishes in 
the same fashion as A — > A c . Does the same thing happen 
in the present case? This is by no means obvious in 
Figure [TJJ which shows the ratio Si p /S dropping slowly 
as A increases, but nowhere near zero. 
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TABLE II: Series coefficients of A^ for the binding energies in the channels S = 0, 1, and antibinding energy (S = 2). The 
S — and S = 2 states are doubly degenerate. 



N 


S = 


S = 1 


S = 1 


S = 2 





0.00000000000000D+00 


0.00000000000000D+00 


0.00000000000000D+00 


0.00000000000000D+00 


1 


1.00000000000000D+00 


5.00000000000000D-01 


5.00000000000000D-01 


5.00000000000000D-00 


2 


-2.25000000000000D+00 


-2.12500000000000D+00 


-3. 12500000000000D+00 


-1.37500000000000D+00 


3 


-1.93750000000000D+00 


1.31250000000000D+00 


-2.93750000000000D+00 


1.87500000000000D-01 


1 


-3.07812500000000D+00 


2.97656250000002D+00 


-2.77343749999998D+00 


2.27343750000000D+00 


5 


3.47656250000001D-01 


1.07812500000003D+00 


3.06250000000002D+00 


2.36718750000000D+00 


6 


-9.69726562500059D-01 


-1.00527343749999D+01 


8.35742187500014D+00 


-8.13476562500000D+00 


7 


3.51385498046887D+00 


7.44207763671879D+00 


4.07301635742189D+01 


-7.26873779296875D+00 


8 




7.92327880859462D+00 


1.69468475341798D+02 


-3.48072814941411D+00 



To pursue this question further, we have studied the 
series at k = (0,0), also given in Table I. A Dlog Pade 
analysis of these series reveals a dominant singularity at 
A = —0.43(1), with exponent around —0.65(3) in both 
cases. This will correspond to another critical point of 
the model, where the spins order ferromagnetically in the 
planes, and antiferromagnetically between them. At pos- 
itive A, there is no sign of a pole around A = 0.4. The 
ratio Sip/S decreases smoothly to around 0.80 at the 
critical coupling, and shows no sign of vanishing there. 
Thus it appears that in this case the renormalized residue 
function does not vanish at A c , except at the Neel mo- 
mentum. 



B. Two- particle excitations 

We have generalized the computer codes which were 
previously used to calculate 2-particle perturbation series 
for 1-dimensional models [3j to cover the two-dimensional 
case. Figure QT] shows the dispersion diagram estimated 
from the perturbation series for 2-particle states at A = 
0.1. We have zoomed in on the region where the bound 
states occur. It can be seen that S = singlet and S 
= 1 triplet bound states emerge below the continuum 
near k = (tt,tt), and S = 2 quintuplet antibound states 
appear above the continuum, as predicted by the triplet- 
wave theory. The S = and S = 2 states are doubly 
degenerate at k = (ir, ir). All states merge with the con- 
tinuum at some finite momentum point k, and for the 
most part they appear to merge at a tangent, as in the 
one-dimensional case [34[ • The results look very similar 
to the triplet- wave predictions shown in Figure [S] 

Figure [T2] shows the behaviour of the binding energies 
at k = (ir, it) as functions of A, as estimated from Pade 
approximants to the series given in Table [ill The degen- 
erate pair of singlet bound states are the lowest over most 
of the range, but merge back into the continuum some- 
what before the critical point. One of the triplet states 
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FIG. 11: Series estimates of the energies of 2-particle states 
at fixed A = 0.1, along symmetry lines in the Brillouin zone. 



disappears into the continuum quite early, but the other 
appears to disappear only at the critical point. For the 
AHC, the binding energies also vanished at the critical 
endpoint of the dimerized phase. The pair of antibound 
quintuplet states, on the other hand, appear to remain 
above the continuum even at the critical point, from our 
estimates. 
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0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 



X 

FIG. 12: Binding energies as functions of A, relative to the 2- 
particle continuum. For bound states, we graph {E2 P — E~ ont ), 
while for antibound states we graph (E2 P —E^ ont ), where E^r ont 
denote the upper/lower bounds of the continuum. 

IV. SUMMARY AND CONCLUSIONS 

In this paper, we have used a modified triplet-wave 
theory and dimer series expansions to study the Heisen- 
berg bilayer system in the dimerized phase. As found in 
earlier papers [13, [IH, [T3|, the model displays a quan- 
tum phase transition from the dimerized phase to a Neel 
phase at a coupling ratio A c = 0.394(1), with critical in- 
dices in good agreement with the predicted values from 
the classical 0(3) nonlinear sigma model in three dimen- 
sions, v = 0.707 and -q = 0.036. 

Our modified triplet-wave approach is found to give 
good results at small couplings A, but towards the crit- 
ical region the self-consistent Born approximation ap- 
proach of Kotov et al. 0, which includes some 
important higher-order terms, gives much better results. 
The triplet- wave approach predicts, as for other dimer- 
ized systems, two-particle bound states in the 5 = and 
5 = 1 channels where an antiferromagnetic alignment of 
spins can give rise to an attractive force, and antibound 
states in the 5 = 2 channel, where the spin alignment is 



necessarily ferromagnetic and repulsive. 

Our series calculations focused upon two major fea- 
tures, the critical behaviour of the static transverse struc- 
ture factor, and the spectrum of 2-particle bound states 
in the model. The integrated structure factor 5(k) and 
the single-particle component 5i p (k) were both found to 
diverge at the critical point for momentum k = (7r,7r), 
with exponents agreeing well with the predicted value 
(rj — 1)^ = —0.68. The ratio Si p /S remains finite 
throughout the region, even at the critical coupling A c . 
This is in contrast to the case of the alternating Heisen- 
berg chain, where the 1-particle component vanishes log- 
arithmically at the critical point [3l|,[35j]. In fact, here the 
one-particle state dominates everywhere (Si p /S > 80%). 

In the 2-particle sector, a pair of bound states is found 
in the 5 = and 5 = 1 channels near momentum k = 
(tt, 7r), as predicted, and a pair of antibound states in 
the 5 = 2 channel, the pairing being a two-dimensional 
effect. The singlet 5 = states have the lowest energies 
at small couplings, but both 5 = states and one 5 — 
1 state merge back into the continuum as A increases, 
leaving only one remaining triplet bound state, which 
appears to merge with the continuum right at A = A c . In 
the S=2 channel, both antibound states appear to remain 
above the 2-particle continnum at all couplings A > 0. 

As one moves away from k = (tt,tt), the 
bound/antibound states eventually merge into the con- 
tinuum also. They appear to merge with the continuum 
at a tangent, much as in the one-dimensional case [3l| . 

In future work, we hope to perform similar calculations 
for other two-dimensional models, such as the simple 
Heisenberg model on the square lattice, and the Shastry- 
Sutherland model, which has already been studied by 
Knetter et al. 0], and where the two-particle states dis- 
play some intriguing behaviour. 



Acknowledgments 

This work forms part of a research project supported 
by a grant from the Australian Research Council. We 
are grateful to the Australian Partnership for Advanced 
Computing (APAC) and to the Australian Centre for Ad- 
vanced Computing and Communications (ac3) for com- 
putational support. 



APPENDIX A 

The vertex functions $4 are: 

$ 4 1) (1234) = -[(71 + 7 2 )(cic 2 + cis 2 + sic 2 + s 1 s 2 ){c 3 s 4 + s 3 c 4 ) + (73 + 7 4 )(sic 2 + c 1 s 2 )(c 3 c 4 + c 3 s 4 + s 3 c 4 + s 3 s 4 ) 
+7i+3(cis 3 - sic 3 )(c 2 s 4 - s 2 c 4 ) + 7i +4 (cis 4 - sic 4 )(c 2 s 3 - s 2 c 3 )] (62) 

$ 4 ; (1234) = -[(71 + 7 2 )(s 3 c 4 + c 3 s 4 )(cic 2 + cis 2 + sic 2 + sis 2 ) + (73 +74)(cis 2 + sic 2 )(c 3 c 4 + c 3 s 4 + s 3 c 4 + s 3 s 4 ) 
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-7i-4(cic 4 - sis 4 )(c 2 c 3 - s 2 s 3 ) + 7i_ 3 (cic 3 - s 1 s 3 )(c 2 c 4 - s 2 s 4 )] 



(63) 



$i 3) (1234) 



(71 + 73)(C1C 3 + C1S3 + S1C3 + SlS 3 )(c 2 C 4 + S 2 S 4 ) + (72 + 74)(C2C4 + C 2 S 4 + S2C4 + S 2 S4)(ciC 3 + Si S3) 
+7l-4(ciC 4 - S1S4)(S2S 3 - c 2 c 3 ) 4- 7i+2(cis 2 - sic 2 )(s 3 c 4 - C3S4) (64) 



$i 4) (1234) 



(71 + 72)(C1C 2 + C X S 2 + SlC 2 + SlS 2 )(c 3 C 4 + S 3 S 4 ) + (73 + 74)(C1S 2 + SiC 2 )(c 3 C 4 + C 3 S 4 + S3C4 + S 3 S 4 ) 

+7i-4(cic 4 - sis 4 )(c 2 s 3 - s 2 c 3 ) + 7i+3(c 2 c 4 - s 2 s 4 )(cis 3 - s x c 3 ) (65) 



We have 'symmetrized' these expressions with respect to their indices, using momentum conservation. 



APPENDIX B. Scaling Theory for Structure 
Factors 

Let us briefly review scaling theory in the vicinity of 
a quantum critical point for quantum spin models on a 
lattice. Firstly, the integrated or static structure factor 

urn 



(66) 



is just the Fourier transform of the spin correlation func- 
tion in the ground state, where S" represents the a com- 
ponent of the spin operator at site j. In the continuum 
approximation near the critical point, this reduces to 

5 Q/3 (k) = / d n re ik r < S a (r)S p (0) > (67) 



where n is the number of spatial dimensions. 

The oscillating factor exp(ik • r) will kill off the contri- 
butions from large distances unless it is compensated by a 
corresponding oscillation exp(— iko • r) in the correlation 
function. Then we can write 



5 Q ^(k) = / d n r e ici r g(r) 



(68) 



where q = k — ko, and g(r) is a smooth function. Scaling 
theory [381 ] then tells us that in the vicinity of the critical 
point 



-(d-2+ v ) 



fir/0 



(69) 



where d = n + 1 is the number of space-time dimensions, 
and £ is the correlation length. Thus when k = ko, the 
'critical momentum', we have 



S Q/3 (k ) = J d d -h 



C 1 -" / d d -h z-^- 2+, ^f(z) (70) 



where z = r/£. As the coupling A — > A c , corresponding 
to a quantum phase transition, we expect 



^ ~ (A c - A)- 



(71) 



and hence 



(72) 



as noted in the text. 

For q small but non-zero, |q| <C l/£, we have 



S a/3 (k) ~ £ 



d a ~ L z z -^- 2 +^ e l ^ ms ^f{z) 



> / d^z' z '-(d-2+,) e «'cos(e)yr/^/) 



(73) 

so that at the critical coupling we expect 5 a ^(k) to scale 
like at small q. 

For the 1-particle structure factor, we may paraphrase 
Sachdev's argument [36| as follows. Assuming relativistic 
invariance of the effective field theory, which applies to 
many though not all models, the dynamic susceptibility 
in the vicinity of a quasiparticle pole is expected to have 
the form 



x(k,w) = 



.4 



c 2 k 2 



(74) 



where e is a positive infinitesimal, c the quasiparticle ve- 
locity, A is the quasiparticle energy gap, and A is the 
"quasiparticle residue" . Then the dynamic structure fac- 
tor is 



1 



/m{x(k,w} 



Let 



E(k) = \/c 2 k 2 + A 2 



(75) 



(76) 



then from (jT4"|) , ([75]) and (|76[) we can write the dynamic 
structure factor for the 1-particle state 



Si p (k,«) 



A(k) 



2E(k) 

and hence the static structure factor 



S(uj - E{k)) 



5i P (k) 



du>Si p (k, oj) 



2E{k) 



(77) 



(78) 
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where A(k) is the residue function. 

From renormalization group theory (38j . the sca ling 
dimensions of these quantities are expected to be [33| 
dim[x] = — 2 + r\ and dim [A] = 77, or in other words we 
expect near the critical point 

A(ko) ~ (A c - X) vu , 

E(k ) ~ (A c -A)", (79) 



and hence 

5i p (k ) ~ (Ac - A)-^", (80) 

just as for the total structure factor. This is the result 
quoted in the text. 
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